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A fast Poisson solver for unsteady incompressible 
Navier-Stokes Equations on the half-staggered grid 


Abstract. In this paper, a fast Poisson solver for unsteady, incompressible Navier-Stokes equations 
with finite difference methods on the non-uniform, half-staggered grid is presented. To achieve this, new 
algorithms for diagonalizing a semi-definite pair are developed. Our fast solver can also be extended 
to the three dimensional case. The motivation and related issues in using this second kind of staggered 
grid are also discussed. Numerical testing has indicated the effectiveness of this algorithm. 
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1. Introduction. Consider the unsteady incompressible Navier-Stokes equations 
(INSE) 


( 1 ) 

(2) 


<9w , Jw i ,,<9w 

-w + u ^ + v ~ 


-f grad p = a div grad w 
div w = 0 


where w = (u,v)', in a two-dimensional region fl with initial condition 

w(i,|/,0) = w°(i,t/) in fi 

satisfying the constraint condition (2), and boundary condition 

w(x,y,t) = v/ b (x,y,t) on dVl. 

The boundary condition on <9f 1 satisfies the consistency condition 

(3) <£ w n ds — 0. 

Jdn 

The main difficulty in the solution of this problem is that the unsteady INSE is not 
entirely evolutionary; it is subject to the divergence-free constraint (2). The projection 
method, proposed by Chorin [3] and Temam [20], proved to be very effective and its 
variants are widely used in many applications. With this method, some auxiliary ve- 
locity is found and then projected into the divergence-free space via the solution of a 
Poisson equation for some form of the pressure. For example, let the discretization of 
(1) and (2), after linearization of nonlinear terms, be 

Aw 


(4) 

(5) 


At 


+ /l Aw + V</> — E 
V • w = 0 


where Aw = w n+1 — w n , A represents the discretization of convection and diffusion, 
V and V- the discretization of grad and div respectively, <f> — p n+1 — p n the pressure 
increment, and E includes all known terms at the nth level. With an approximate 
factorization, Eq. (4) is put into the following fractional steps 


( 6 ) 


Aw 

*A t 


+ y4Aw = E 


where Aw = w — w n (w is called the auxiliary velocity) and 


(7) 


w n+1 — w 


At 


+ V<f> = 0. 


Now, applying the discrete divergence to (7), we obtain the discrete Poisson equation 
for <j): 

( 8 ) 

in which the discrete divergence free condition (5), V • w n+1 = 0 has been enforced. 
The solution process is then: 
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Fig. 2. Half-Staggered Grid 


• Find w by (6). 

• Solve the discrete Poisson equation (8) for <j>. 

• Update w n+1 by (7). 

For a more detailed account with the boundary condition treatment , see [10]. 

The constraint (5) causes difficulties in the choice of the mesh, the coordinate 
system and the related discrete Poisson equation, etc., of the solution method. For a 
rectangular region fl, the usual staggered mesh (u, v and p staggered mesh as shown 
in Fig. 1) is most frequently used. With this mesh, no pressure boundary condition is 
needed, which is mathematically correct. Let the discretized linear system of equations 
of (8) be denoted as 

(9) = R , 

where $ is the unknown vector of size, say, M. With this mesh, L is of rank M — 
1; the constraint on the right hand side R reduces to an approximation of (3), and 
the solution $ is unique up to a constant. For a rectangular region with uniform 
or nonuniform mesh intervals, direct solvers from the FISHPACK have proved to be 
most efficient, see [2, 16] for example. Also for the usual staggered mesh, the finite 



difference approximation of (1) is the same as the finite volume approximation. In 
addition, Vp and V • w are straightforward without any interpolation. However, u and 
v involve different finite volumes and near the boundary half interval differencing is 
necessary, which are inconveniences. More serious is the problem on curvilinear meshes 
in arbitrary regions; there the extension of the usual staggered mesh will lose many 
of the above advantages, see [17] for instance. Hence, the half-staggered mesh (w, p 
staggered as shown in Fig. 2 was studied in [10, 11]. This mesh was first proposed 
in [5] for the unsteady INSE; it was recently used in [1] with a different discretization 
of projection, for example; but it is not used in practical computation in general. For 
steady-state INSE, this mesh was used and investigated for FEM in [15] and for the 
Galerkin formulation in [18]. 

The advantages of this mesh for the unsteady INSE are: no pressure boundary 
condition is required; the finite difference approximation of (1) is the same as the finite 
volume method with the same finite volume for u and v; Vp and V • w involve only simple 
averages, and there is no half interval differencing near the boundary. Most important 
is the extension of this staggering to curvilinear meshes in arbitrary regions. As the 
velocity components are at the same point, it is easily transformed to, say, contravariant 
components for finite volume formation, for approximation of div w, and for boundary 
conditions, etc.. The disadvantage of this mesh is the nature of the discrete Poisson 
equation (8), where L is of rank M — 2. There is an additional constraint and the 
solution $ may have oscillations; see the next section. 

As the first step to this problem, we consider rectangular regions with nonuniform 
meshes. Since (9) is to be solved at every time step for the unsteady INSE, efficient 
solvers are necessary and this is the subject of this paper. 


2. Discretization. Consider a rectangular region with a nonuniform mesh, schemat- 
ically show in Fig. 3. We generate the nonuniform mesh by some smooth transformation 
x(C), p(p), with a uniform mesh in the computational (tj region, and approximate 

d_ _ d(d_ 6_ _ AC 

dx dx dC, ^ Sx 6x AC 


where 8 on the right hand side denotes centered differencing, and similarly for d/dy. 
Then for every interior w point 


( 10 ) 


v ^t+l.*+4 = 


cU 1 

l 

m 

. w . 

“ 2 


^j+i.t+i — . ^2±lA — !ti£ 

X 3 + 1 X j X j + 1 X j 

^j + l.fc+l ~ 4>i+\,k . <^J,A :+l ~ <7 ^, k 
Vk+i - yk ~ r Vk+\ ~ Xk . 


and (7) on all the interior w point can be written as 


( 11 ) 


W n+1 - W 
At 


+ = 0 


where G is the gradient matrix for the solution vector 4>. For every interior (f> point, or 
rather every ( j , k ) cell, 


1 


? 


i 


2 x-, \ — x,_i 
2 J 2 






Fig. 3. Non uniform mesh and the stencil 


2 ?/fc+i 




which results in a matrix form 


(13) DW = -B. 

Here D is a matrix form of V- with boundary modifications, i.e. for cells adjacent to 
the boundary, some of the elements in w in (13) are known and put into B. Hence, the 
final discrete system of equations is: 

(14) L<D = DG$ = R = (DW + B) 

After some tedious derivation, the matrix L can written as 


L = DG = [C n ® D x A m + D y B n ® C m ] , 

which is a nine-diagonal unsymmetric matrix. Let dx{ = x, +1 — x,, dy, = y, +1 — j/j 
and dx, t , + i = (dx,+i + dx,)/2 and dyj il+ j = (dy, +1 + dy,)/2. Then the matrices in the 
expression of L can be written as 


D x 


Dy 


1 

dx\ 

dx 2 


1 

^2 



1 


nXn 





1 1 
1 2 1 


Ci = 


2 1 
1 1 


I Ixl 


j/ l T ) ! 


1 

c/xj2 


1 - f 1 + M 

dzi 2 V^x 12 T da: 23 / 


1 

d*23 


dx. 


_1 / 1 , 1 

i — 2 , m — t \ d-T 7n _ 2,m'— 1 d'X 7n — l i7 




J 

m — 1 ,m 


l 

dx Tf j _ i 
1 

dx-rn— 1 


B„ = 


1 

rfyi2 


l 

dj/12 


1 _ / 1 . l 'l 

^2/12 V^12 ^2/23 / 


1 

3^23 




_I / 1 , 1 

Vi— 2,n— 1 \ ^Vn— 2,n— 1 dj/ n _ 1 ,n 


dj/n— l,n 


1 

dj/n- l,n 

1 

dyn — l t n 


Here an m x n rectangular grid is assumed (Note that the n used here is different from 
that used in the time stepping process.) For the interior nodes, a nine-point stencil is 
presented in Fig. 3. The corresponding coefficients for a node ( z , j) in that figure are: 
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11 11 
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d{J — 


(fxj i dyj dyj—ij 
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dx{ \dx,-i'i dx 

11 11 
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1 \ 1 1 

m+i / dyj dyj-ij 


dx{ dxi i^-i dyj dyj—\ j 

n 1 1 1/1 1 

2 I 1 
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11 11 
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l \ +2 J 1_ 

.,.+1 / dyj dy h j + i 



1 1 
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°iyj “ 
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1 1 


dx± dx dxjj dyj j^.\ 
( a t,i “f“ C i,j d“ fi,3 “t" ^t,j) 


It is not difficult to see that L is singular and in the next section we will show it is 
of rank M — 2. The two base vectors of the null space N(L) are 

4*oi =(1, 1,1* *, i) T ®( 1,1, 1---,1) T 

— e n 

4*02 = (1, — 1, 1 * • • , — 1") T <S> ( — 1, 1, — 1 * ’ * , — l m-1 ) T 

— G n ® ® m 

where 

e, =(1,1,1,...,1) T 

S v- ' 

l 

e'i =(l,-l,l,---,-l'- 1 ) T 
' ' 

l 

since A m e m = 0, B n e n - 0 and C/e'; = 0. 

We note from 4>oi that the solution 4* can differ by a constant which is irrelevant 
as our interest is only in V^. From 4> 02 , we see that 4> can have oscillations, called the 
checkerboard effect [15], and from (10), we also see that in the present finite difference 
method the oscillations in 4> do not affect V 4> and hence the solution w. 

To ensure the solution of (14) exists, the right hand side R must satisfy 

R e N{L t ) ± . 

Thus, R has to be orthogonal with the two base vectors of N(L T ). For a nonuniform 
mesh, the two base vectors of the null space of the discrete operator ZT are 

V>oi = dy n <g) dx m 
^02 = e' n ® —e' m 

where 

— [dx\j dx 2 , * * * 1 dx m ) . 

V v ' 

m 

and 


dy„ = (dy u dy 2 ,- ■ ■ ,dy n ) T . 

S . . > 


71 



This can be verified by multiplying directly with the matrix L T . 

L T ip 01 = [Cn ® Am D x + B n Dy ® Cm] dy n ® dx m 

= C n dy n ® A m D x dx m + B n Z) y dy n ® C m dx m 
= 0 

L T xp 02 — [Cn ® A m D x + B n D y ® ( 7 m ] e' n ® — e' m 

C n & n m B n DyG n ® ^ m 

= 0 

since y4 m Z) x dx m = 0, B n D y dy n = 0 and C/e'/ = 0. 

In [15], Sani investigated the solution of (13) and gave the base vectors of N(D T ) 
to be xpoi and i/> 02 - He showed that the constraint ip^B — 0 reduces to a discrete form 
of (3), and ip$ 2 B = 0 yields a constraint for the tangential velocity component on the 
boundary. 

For the unsteady INSE, the interest is in the solution of (14). Due to the divergence 
form of /?, ip^ x DW will cancel out, leaving ip^ R = ( xpQ 1 B)/At . Hence, the constraint 
R = 0 is a discrete form of (3). ipQ 2 R = 0 will a ^ so hold if (5) holds for some velocity 
field with the given boundary condition, i.e. (13) holds for some VF, for example, the 
initial velocity field with time independent boundary condition. 

When dxi = Ax and dy t = Ay are constant, the matrix L is symmetric and has a 
simpler form 

i= + 

where 

' -1 1 
1 -2 1 

4 = 

1 -2 1 
1 - 1 

In particular, when Ax, = Ay, = h, the discretization is a well known skewed five point 
stencil (see Fig. 4). 

3. Fast solver. As we mentioned in the introduction, since the matrix equation 
(9) must be solved every time step, an efficient solver is essential. Fast solvers for 
the sum of two matrix tensor products were discussed in the first author’s paper [2], 
Later, L. Kaufman and D. Warner [13] provided a generalization to the higher order 
discretization case. In particular, the eigen-decomposition applused in [2] is replaced 
by a generalized eigenvalue and eigenvector problem. If we examine the formula of the 
matrix L, it is clear that generalized eigen-decomposition is needed to develop our fast 
solver. For a matrix of the form 



A ® B + C ® E, 



2 2 

• o • 

• o • 

2 2 

Fig. 4. Skewed five point stencil 

the key step in developing a fast algorithm in [2, 13] was to simultaneously diagonalize 
the matrices B and E by a matrix Z for which 

Z t BZ = I, Z t EZ = D, 

where B and E are symmetric and positive definite. There are two obvious problems 
which cause these algorithms not to be directly applicable here. The first problem is 
that the matrices L , D x A m and D y B n are not symmetric. The second problem is the 
singularity of the matrices A m , Ci and B n . Many of the effective generalized eigenvalue 
algorithms [4, 8, 12, 21] for diagonalizing two banded matrices simultaneously need 
the matrix pair to be nonsingular and symmetric. Unfortunately, both A m and C\ are 
singular and D x A m is unsymmetric and thus a new approach is needed. 

3.1. Equal-spaced in one direction. Assume the mesh is equal-spaced in one 
direction, say the x direction. 

Let A X{ = Ax , i = 1, • • • ,m. Then the matrix L can be written as 
(15) L = — — 2 C n (g) A' m + D y B n ® C m . 

Though the matrix L is not symmetric, both — A' m and C m are symmetric, semi-definite 
matrices. We can show the following result. 

THEOREM 3.1. There exists a matrix Z such that 

Z T C 2 Z = —A' m 
Z T S 2 Z = C m 

where the diagonal elements Ci, s, of C and S satisfy 

c i + s i = 1> 0 5: c i, -s, < 1 i = 1 , • • • m. 

Proof. The proof is a natural generalization of the result presented in [21], where a 
pair of symmetric and definite matrices is considered. 



There exist two orthogonal matrices T x and T 2 such that 

-A' m = T X D\T X , 

C m = T 2 D 2 2 Tj, 

where D,, i — 1,2 are diagonal matrices with non-negative diagonal elements. Now, 
construct 

W — _ 

w “ [ d 2 tj J “ 

where Q , R make up the QR decomposition of the matrix W. Then partition matrix 
Q as follows: 

(. 6 ) «=[£’ 

Applying Stewart’s theorem[19], we have from the Singular Value Decomposition(SVD)[7] 

Qx = U X CV X 
Q 2 = u 2 sv 2 

where U x , U 2 , V\ and V 2 are all orthogonal matrices and C, S are diagonal matrices 
diag(cj) and diag(s t ), respectively. Then from (16) and the orthogonality of Q, 

v x = v 2 = v 

and 


Therefore, we have 

-A' m = T X D\T X = R T V T C 2 VR = Z T C 2 Z, 
c m = t 2 d\t. J = R t V t S 2 VR = Z T S 2 Z , 

where 


Z = VR. 


□ 

After the simultaneous diagonalization of the matrices A m and C m is achieved, 
the fast solver in this case is similar to the Buzbee-Golub-Nilson or Kaufman- Warner 
algorithm. First, the matrix equation (15) can be transformed into a tridiagonal matrix 
equation as follows: 


(I®Z~ T ) L (/0Z- 1 ) 


/®Z- ~T^ C n ® A' m + DyB n ® C m I®Z~ X 

-^C n ® Z~ T A' m Z~ x + D y B n ® Z~ T C m Z~ x 

c n ® C 2 + D y B n ® S 2 

L", 



where L" is a tridiagonal matrix. A simple permutation matrix transforms L" into a 
block diagonal matrix. The diagonal elements are n x n tridiagonal matrices. The 
solution of these tridiagonal matrices is trivial. In particular, they can be implemented 
on a vector or parallel computer efficiently. 

For a regular mesh, Ax = Ay = h , this block diagonal matrix has a very simple 
form. We show the following result. 

THEOREM 3.2. The rank of the matrix L for regular mesh is of mn — 2. 

Proof. We know that both A' m and C m have only one zero eigenvalue and N{A' m ) fl 
N{C m ) = 0. Therefore, the two zero diagonal elements c,. and Sj. in C and S will not 
be in the same location. It is easy to arrange them such that cj = 0 and s m = 0. After 
permuting the matrix L" to a block diagonal matrix, the first and last diagonal block 
matrix will be: 

T, = S 2 A' T — r 2 r 

J n — C n <_/ n , 

respectively. Each has only one zero eigenvalue. The rest of the diagonal block matrices 
are of the form 

( 17) Ti = —c 2 C n + s 2 A' n , i — 2 , • • • n — 1 , 

where c 2 + s 2 = 1. When c, = s t , (17) reduces to a diagonal matrix 

diag{ 1, 2, • • • , 2, 1}, 

and it is not singular. If c, ^ s,, Ti is a tridiagonal matrix, viz 

—a 1 
1 —2a 1 

T t = {-c 2 + s 2 ) 

1 -2a 1 

1 —a 

where |a| = — j— i — j- > 1, and hence T, is non-singular. □ 
c t 

We conjecture that this is true in general. Numerical computations support our 
conjecture. 

As we discussed in the previous section, the solutions of the matrix equation (14) 
can differ by a linear combination of the two base vectors in N(L). Some special 
attention needs to be paid to the solution of linear systems involving the first and last 
tridiagonal submatrices 7\ and T m . 

3.2. General case. For a general non-uniform grid, the simultaneous diagonal- 
ization of both matrix D x A m and C m by a single matrix can not be achieved. However, 
the generalized eigenvectors 1 can still be used to achieve the same purpose with two 

1 It is not difficult to see that the generalized eigenvalue problem: 


&xA m x — XC m x 



different matrices. Recall the matrix equation: 

(19) L$ = R 
where 

(20) L = DG = [C n 0 £> r A m + £> y 5 n ® C m ] . 

The QZ algorithm [7, 14] is used to compute the generalized eigenvalues (cq,/?,) and 
its corresponding eigenvectors v, such that 

0i D x -4mW. — ^iCrn^t* 

Let T> q and "D c be two diagonal matrices whose diagonal elements are a, and /?,, respec- 
tively. We have 

D x A m VV c = C m V2> a 

where matrix V = [ui, Ui, • • • , u m ]. Since we know that each matrix A m and C m has one 
zero eigenvalue, it is clear that each of T> a and T) c has one zero diagonal element. Due 
to the orthogonality of the null spaces of N(A m ) and N(C m ), the two zero diagonal 
elements will not be in the same location. We arrange the generalized eigenvectors such 
that 


i = 0; A = 0. 


Then we construct a matrix 

(21 ) U — \D x A m [i>i , t>2 , ■ ■ ■ .u m _ i], C m %] . 

It is not difficult to see that 

' 1 


and 


U~ 1 (D x A m )V = 


1 


= L 


02/ q 2 


U-'CV = 


= V X . 


fin-i/a. 


n— 1 


is a discrete approximation of the eigenvalue problem 

(18) - y" = Ay, y'(0) = y'(l) = 0 

on a non-uniform grid. Therefore, this pencil is diagonalizable. Our extensive numerical testings with 
randomly non-uniform grid have verified this observation. The eigenvectors are approximations of the 
eigenfunctions of ( 18 ). 



Thus, the tridiagonalization of (20) can be achieved as follows: 


(JnOtr 1 )^,® V) = {In®U-')[C n ®D x A m + D y B n ®C m ){I n ®V) 

= C n ®U~'D x A m V + D y B n ®U~ l CV 
= C n ® 7 m + D y B n ® V x 
= L. 

Here, the matrix L has three non-zero diagonals with a bandwidth of m. The resulting 
matrix equation will be 

(In®U- 1 )L(I n ®V)(I n ®V- 1 )$ = ( I n ®U~ l )R 

L$ = R , 

where 

$ = (/^y- 1 )^ 

R = (In^U-^R. 

We start with an ordering of the grid nodes first in the x direction, then in the y 
direction. Let P be the permutation matrix, which reorders the grid nodes first in the 
y direction, then in the x direction. We have 

(. PLP)(P * ) = PR 
= R, 

where L = PLP , $ = P4> and R = PR. The matrix L is a block diagonal matrix in 
which each diagonal element T, is an n x n tridiagonal matrix: 

T i = C n 

T, = C n + — D y B n (* = 2,3, • • ■ ,m — 1) 

CX{ 

T m = D y B n 

and 



The solutions of the T/s are trivial, but, some attention should be paid to the first 
one and last one, since both T x and T m are rank n — 1 and the top left (n — 1) x 
(n — 1 ) submatrices of T\ and T m are not singular. We can solve the two non-singular 
submatrices and assume the last element of the solution is zero. Using this approach, 



the contribution of the two null space base vectors of L is set to be zero and this will 
not affect the final solution of the INSE. The final solution of (19) is 

$ = (/„ ® V)P$. 

The fast solver can be summarized as follows: 

1. Preprocess. 

(a) Compute the generalized eigen-decomposition {T> a ,'D c ,V) of the pencil 
( D x A m , C m ) • 

(b) Compute the LU decomposition of the matrix U given in (21). 

(c) Form the LU decompositions of the matrices T{. 

Note that this step is needed once at the beginning of the solution process of the 
INSE. The cost of it can be amortized over many time steps. The complexity 
of step (a) and (b) is 0(m 3 ), and of (c) is 0(nm). 

2. Solve 

(/ ® U)R = R. 

The complexity is 0{nm 2 ). 

3. Solve 

L$ = R = PR. 

The complexity of this step is O(nm). 

4. Compute 

$ = (I ® V)P$ 

The complexity of this step is 0(nm 2 ). 

4. Three dimensional case. The fast solver can be extended to three dimen- 
sional problems without difficulty. Assume an m x n x l non-uniform grid for the 
solution domain. The coefficient matrix of the three dimensional case in (19) will be 

L = Ci® C n ® D x A m + Ci® D y B n ® C m D Z E{ ® C n ® C m . 

where 

1 

(Til 

(Tz^ 

1 

(Ti/ J /x/ 




and 


dz\ 2 dz\ 2 

d^i2 ( d^i2 


There are two generalized eigenvalue problems 

BxA m X — \C m X , 

D y B n y = 7~C n y 

which need to be solved. Let V and IT be the eigenvector matrices for the two problems, 
respectively. Similarly, form 

U [D x A m [v \ , * * * •U 7 n _i], CmVjji], 

and we have 

U~ l D x A m V = I m , U~'C m V = V x , 

= 7 n , (5- 1 C„VF = Py. 

The tridiagonalization of T can be achieved as follows 

(/( 0 <? _1 0 Im)(Il telnttU- 1 ) L (I,®W® I m )(h ®I n ®V) = 

Cl 0 Vy 0 7 m + C, 0 /„ 0 V x + D Z E, ®V y ®V x 

5. Numerical Experiment. Consider INSE with exact solution 

u — e t sin x cos y 

(22) v = — c* cos a: sin y 

p = e l sin x sin y 

with corresponding nonhomogeneous terms in the momentum equations, on a square 0 < x < 
7 r and 0 < y < x. The initial and boundary values of u and v are taken from (22). The 
conservative version of the centered-difference Crank-Nicolson scheme was used for (6). For 
the discrete Poisson equation, the first constraint 4>q X R = 0 is obviously satisfied, allowing for 
discretization error. The second constraint <$ 2 R = 0 reduces to 

(23) £v w -£Vw = 0 

red black 

corresponding to the checkerboard effect, and hence is also satisfied (except for the discretiza- 
tion error). 




On a 64 x 64 nonuniform mesh with mmi(dxi) — min i(dyi) « 3.3 X 10” 2 and ma X{(dxi) — 
ma Xi(dyi) as 7.6 X 10~ 2 , the fast solver developed herein took 0.13 sec CPU time on the 
100MHz Indigo R4000 workstation for each solution of (14), with residual % 0.4 • 10 -3 and 
max t j(V • w) % 0.1 * 10 -3 . There was no apparent oscillation in pressure, but adding ±1, 
say, to the initial pressure at alternate cells (the red and black cells) will cause an oscillating 
pressure distribution as shown in Fig 6 for t = 1. As we indicated earlier, the oscillation in p 
has no effect on the solution of u and v shown in Fig. 5, nor any effect on the solution process. 

In contrast, solving the same problem using an adaptive multigrid method [10] took six 
times more CPU time. This version of the multigrid is not very effective in this application, 
presumably due mainly to the treatment for the second constraint. 



Fig. 5. Solution u and v at t — 1 

6. Conclusion. We have developed a fast Poisson solver for the unsteady incompressible 
Navier-Stokes Equation with finite difference methods on the non-uniform, half-staggered 
grid. Due to the efficiency of this method, it is possible to solve the unsteady flow with 
Re = 10,000 [9]. Although this method can only be applied to an orthogonal rectangular 
grid, it can be used as a preconditioner or in a domain decomposition scheme for general 
applications. Our next project will be development of effective solvers for unsteady INSE on 
an irregular domain with curvilinear grid. 
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